The fate of phonons in freely expanding Bose-Einstein condensates 
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Phonon-like excitations can be imprinted into a trapped Bose-Einstein condensate of cold atoms using light 
scattering. If the condensate is suddenly let to freely expand, the initial phonons lose their collective character 
by transferring their energy and momentum to the motion of individual atoms. The basic mechanisms of this 
evaporation process are investigated by using the Gross-Pitaevskii theory and dynamically rescaled Bogoliubov 
equations. Different regimes of evaporation are shown to occur depending on the phonon wavelength. Dis- 
tinctive signatures of the evaporated phonons are visible in the density distribution of the expanded gas, thus 
providing a new type of spectroscopy of Bogoliubov excitations. 



PACS numbers: 03.75 Fi 
i I. INTRODUCTION 

Bose-Einstein condensates of cold atoms represent a re- 
markable testing ground for theories of weakly interacting 
bosons. The fact that the actual condensates are spatially 
confined and inhomogeneous makes these systems even more 
interesting. From one hand, the Bogoliubov theory of uni- 
form gases |1] can still be applied in certain limits, when 
the system can be viewed as "locally uniform". Such a lo- 
cal density approach has been used, for instance, to character- 
ize the response of a condensate to light scattering processes 
. In the same spirit, properties of the uniform 
gas can also be extracted from a trapped condensate by reduc- 
ing the effects of inhomogeneous broadening as suggested in 
Ref. 1 9]. On the other hand, finite size effects are interesting 
by themselves, giving rise to a new type of Bogoliubov excita- 
tions with spatially varying quasiparticle amplitudes, such as 
the low frequency collective oscillations of the whole conden- 
sate (see I IQl and references therein) and the axial phonons 
of e long ated condensates, exhibiting a multibranch spectrum 

In most of the experiments with excited condensates, the 
observations are performed after switching-off the confining 
potential and letting the condensate to expand. The compar- 
ison between theory and experimental data often ignores the 
dynamics of the expansion. This is justified if one looks for 
quantities that are conserved during the expansion like, for in- 
stance, the total momentum of the condensate. However, the 
expansion of a condensate initially dressed with phonon-like 
modes is interesting from both the conceptual and experimen- 
tal viewpoints. Conceptually, the behavior of an excited state 
that starts as a quasiparticle and evolves into a particle is a re- 
markable and nontrivial example of quantum process, which 
shows some similarities with the quantum evaporation pro- 
cess at the surface of superfluid helium and with the evo- 
lution of two-level systems with nonhermitian Hamiltonian 
fTiil . From the experimental viewpoint, on the other hand, the 
characterization of the observable density and velocity distri- 



butions of the expanded gas in terms of specific initial config- 
urations is important in order to use those distributions as a 
probe of in-trap quasiparticles. 

In this paper we use the Gross-Pitaevskii (GP) theory to 
investigate the expansion of elongated axially symmetric con- 
densates. We assume the initial configuration to be a station- 
ary trapped condensate at zero temperature. The condensate 
is excited by populating some quasiparticle states and, then, 
let to expand. 

In section II, we first present the results of numerical sim- 
ulations based on the direct integration of the time dependent 
Gross-Pitaevskii equation in the case of an elongated con- 
densate similar to the one of recent experiments |(4|, [Till . 
We show that the expansion exhibits quite different behav- 
iors depending on the wavelength of the initial phonons: i) 
long wavelength phonons remain inside the expanding con- 
densate in the form of density modulations; ii) short wave- 
length phonons are converted into a separate cloud of excited 
atoms moving out of the condensate. In both cases, the final 
density distribution shows nontrivial features, such as a radial 
distortion of the density modulations in the condensate and a 
"shell" structure of the released-phonon cloud. 

In section III, we use the simplified geometry of an infi- 
nite cylindrical condensate in order to get a deeper insight, 
pointing out the different role of radial and axial degrees of 
freedom. In this case, we numerically solve the Bogoliubov 
equations for the amplitudes of in-trap quasiparticles. Then 
we follow the expansion of a condensate initially dressed with 
one of these quasiparticles by using the scaling properties 
of the GP equation in two-dimensions and solving rescaled 
Bogoliubov-like equations. This approach allows us to char- 
acterize the behavior of the excitations during the expansion, 
at different levels of approximation, and to find the relevant 
timescales for the evaporation process. Interesting results can 
also be obtained by averaging out the slow radial motion of 
the excitations in the rescaled coordinates. The problem is 
thus mapped into the evolution of a quasiparticle in a uniform 
gas with a decreasing time-dependent density. In section IV, 
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we show that this process can be either adiabatic (conversion 
of a quasiparticle into a single particle with the same momen- 
tum) or non-adiabatic depending on the phonon wavelength 
and on the chemical potential of the gas. 

In section V, we compare the predictions of the rescaled Bo- 
goliubov equations for the infinite cylinder with the results of 
GP simulations for elongated condensates. The comparison 
is instructive and allows one to relate the properties of in-trap 
quasiparticles with observable features of the expanded gas. 
In particular, we discuss the axial motion of density modula- 
tions, the radial scaling of nodal lines, the conditions for the 
appearance of a separate released-phonon cloud and for the 
adiabaticity of the quasiparticle evaporation process. 

The expansion of condensates with long wavelength 
phonons was previously studied in Ref. 1 15|, where random 
phase fluctuations were included in the configuration of a very 
elongated condensate at the beginning of the expansion, in or- 
der to simulate the behavior of a quasi-condensate with ther- 
mal excitations. Low energy phonons that expand in the hy- 
drodynamic regime were also considered in ilill as due to 
"quantum vacuum" phase fluctuations. Instead of investigat- 
ing the effects of thermal and/or quantum fluctuations, we are 
here considering situations where the initial excitations are 
imprinted in a controllable way and characterized by looking 
at the shape of the expanded gas, thus allowing for a new type 
of spectroscopic studies of Bogoliubov quasiparticles. 



n. GP NUMERICAL SIMULATIONS FOR ELONGATED 
CONDENSATES 

We simulate a process of quasiparticle creation and expan- 
sion in a realistic axially symmetric, cigar-shaped, conden- 
sate. The starting point is the time dependent GP equation for 
the order parameter y, z,t) of a condensate of N bosonic 
atoms of mass m LIO, 17]: 

thdt-^ = (^-^^ + V + gN\-^\^^^ , (1) 

where g = Aiih^a/m, and a is the s-wave scattering length 
that we assume to be positive. The order parameter is here 
normalized according to J drj^'p = \. 

Let us first take the harmonic trapping potential in the form 

V{p, Z) = V^trap(p, Z) = (l/2)mc^2(^2 ^ ^2^2) (3) 

where p ^ [x^ + y^]^/^ and A = Uz/ujp. The ground state 
can be found as the stationary solution of Eq. Q. In practice 
we map the order parameter on a Np x grid (typically, 
64 X 1024 points) and propagate it in imaginary time with an 
explicit first order algorithm starting from a trial configura- 
tion, as described in Ref. If A < 1 the condensate at 
equilibrium is a prolate ellipsoid. 

Excitations can be created in the condensate by using light 
(Bragg) scattering |2, 3,"5|5llll- Within the GP theory, this is 
accounted for by adding an extra term to the external potential. 
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FIG. 1: Column (below) and linear (above) density, in arbitrary 
units, of a condensate subject to a Bragg pulse with k — O.Sop ^, 
uj = l.ltjp, Vb = Q.2hL0p and tn = 40^^^, and a subsequent free 
expansion oft = 25a;^^. The white ellipse is the shape of the con- 
densate at t = 0, at the beginning of the expansion. Distances are in 
units of ap — [h/{mijjp)]^^^. 




FIG. 2: Same as before, but for a Bragg pulse with k — 2.5ap ^, 
u — 4.7ujp, Vb ~ ITiLOp and t_g = lujp^. 

that becomes iHqIi 

V{p, z, t) = Vtrap(/0, z) + Vb cos{kz - ujt) (3) 

for —tB<t<0 and assuming the momentum transfer to 
be along z. We numerically integrate the GP equation Q by 
propagating the order parameter in real time, starting from the 



ground state configuration at t — —ts- We split the propa- 
gation into the axial and radial parts. The former is obtained 
by using a fast Fourier transform algorithm to treat the ki- 
netic energy term, while the latter is performed with a Crank- 
Nicholson differencing method, as in Refs. |8., 1 1, _12|. 

Now, let's suppose that both the Bragg pulse and the ex- 
ternal trapping potential are switched-off at i = 0. The con- 
densate freely expands and the expansion can still be studied 
by solving the GP equation. We use the same algorithm as 
before with the only difference that we rescale all distances 
during the expansion in order to keep the box size and the 
number of grid points the same as for the trapped condensate. 
For that we use a dynamically rescaled GP equation already 
introduced in I20l 12111 . 

For our simulations we choose the parameters of the con- 
densate used in the experiments of Nir Dayidson and co- 
workers illllilllll, made of iV = 10^ atoms of ^'^Rb in a 
ti-ap with ujp = 27r(220Hz) and A = 0.114. We plot distances 
in units of the harmonic oscillator length ap — [h/ {mujp)]^/^, 
time in units of a;~^ and energies in units of hujp. The chem- 
ical potential is fi = 9.1hujp. Several condensates in current 
experiments have similar parameters. 

Typical results are shown in Figs. [0 and |2l where we plot 
the column density, J dy\'i>{x, y, z)p, and the linear density, 
/ dxdy\'i>{x, ?/, z)p, for two sets of parameters of the Bragg 
pulse. The values are chosen such that the lowest longitudi- 
nal modes are resonantly excited. Moreover, both values of k 
are smaller than ^"^ ~ Sa^ ^, where ^ is the healing length. 
This implies that in both cases the excitations are collective 
phonon-like Bogoliubov quasiparticles. 

There is a dramatic difference between the density distri- 
butions for the two values of k. At low k all atoms remain 
within the expanding condensate that exhibits an evident den- 
sity modulation, while at high k an extra cloud of excited 
atoms separates out in the positive z direction. 

In order to get a better view of the behavior of expanding 
phonons one can repeat each simulation twice for the same 
condensate with and without excitations (Bragg pulse "on" 
and "off) and calculate the density difference An{p, z) = 
N[\^on{x,y,z)\^ - |*off(a;,y,z)p]. When the number of 
quasiparticles is much smaller than N, this quantity converges 
to the density variation, 5n{p, z), associated with the excita- 
tions in the linear response regime. In Figs.|3]and|3we show 
the results for the two sets of parameters of Figs.^andl^J re- 
spectively. 

At fc = 0.5a~^ we observe well defined wavefronts. They 
reflect the density modulations associated with the initial 
phonon att — Q (dashed line in the upper plot). These mod- 
ulations move along z during the expansion. Their velocity 
depends on both z and p and is faster for p = 0, where the 
density is higher. This causes each front to be slightly bent. 
Effects of such a bending are also visible in Fig.[T]in the form 
of a small asymmetry of the modulations of both the column 
and linear densities. In Fig.|5lwe show the position z{t) of 
a wavefront for p — 0, 0.4 and 0.8R (points A, B and C in 
Fig- 01, where R is radial size of the expanding condensate. 
The initial velocity is very close to the sound velocity in a 
cylindrical condensate, c = [gn{0)/{2m)]^/^ L2Z1 . having the 
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FIG. 3: Lower part: the quantity An(p, z) att — 25LUp^ is plotted, 
in arbitrary units, for the same simulation of Fig. Q Black (white) 
means positive (negative) An(p, z). The dashed line is the position 
of the surface of the expanding condensate. The uniform grey colour 
far outside the condensate corresponds to An = 0. Solid lines are the 
position of maxima of An{p, z) (wavefronts). Dot-dashed lines are 
the same wavefronts estimated by means of Eq. <49t . Upper part: in- 
tegrated density variation 2n J dppAn, calculated at t = (dashed), 
10 (dot-dashed) and 25ujp'^ (solid). 




FIG. 4: Lower part: the quantity An{p, z) alt — 25ujp^ is plotted, 
in arbitrary units, for the same simulation of Fig.|2| The meaning of 
the greyscale is the same as in the lower part of Fig.|3| Upper part: 
integrated density variation 2n f dp pAn, calculated al t — 0, 10 
and 25LOp^ from thin to thick solid line, respectively. 



same density profile n{p) at z = 0. 

At fc = 2.5a~^ the main feature is the shape of the cloud of 
excited atoms that are moving out of the condensate. By us- 
ing the order parameter y, z), one can calculate the local 
contribution to the total energy and momentum of the system. 
Most of the energy of the condensate is, of course, in the ki- 
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FIG. 5: Position z{t) of a wavefront for k — 0.5ap ^ and p — 
0, 0.4 and 0.8R (points A, B and C in Fig. 3), where R is the radial 
size of the expanding condensate. The thin solid line is z(t) — ct, 
where c — [gn{0)/{2m)]^^^ is the sound velocity in a cylindrical 
condensate with the same central density n(0). 



netic energy of the fast radial motion. What is more important, 
however, is that the initial excitation energy, i.e., the energy of 
the excited phonons, is almost completely transferred to the 
lateral cloud, which is moving along z at a velocity of the or- 
der of fik/m. The shape of this cloud exhibits an interesting 
"shell" structure with some lobes divided by almost empty re- 
gions. The appearance of this extra cloud is accompanied by a 
strong suppression of density modulations in the condensate. 

By performing several simulations at different k we find 
a transition from the low-fc scenario (density modulations 
within the expanding condensate) and the higher k scenario 
(external cloud of atoms) occurring in between fc = 0.5 and 
la~^. The experiments so far performed with phonons excited 
by Bragg scattering f?, "s", "4", |^ [Tl| belongs to the second 
scenario. Using a tomographic imaging method the authors of 
Ref. 1 5] found that the observed released-phonon cloud has in- 
deed a nontrivial shell-like shape. First evidences of the first 
scenario have also been found in recent experiments of the 
same group |23|, where density modulations have been ob- 
served at low k within the condensate. Before coming back 
to the interpretation of our GP simulations, let us explore the 
instructive case of an infinite cylindrical condensate. 



III. BOGOLIUBOV EXCITATIONS IN EXPANDING 
CYLINDRICAL CONDENSATES 



FIG. 6: Spectrum of axially symmetric Bogoliubov excitations of a 
cylindrical condensate with 77 = 9.1. The frequency ujnk, in units 
of the radial trapping frequency ujp, is plotted as a function of the 
axial wavevector k, in units of a^^ . The number of radial nodes is 
n = 0,1,2,..., starting from the lowest branch. The quasiparticle 
amplitudes u and v of the modes a, b, c, d, e and / are shown in the 
first column of Figs.Qand[8| 



and obeys the stationary GP equation 



._JL + Vip)+gnoip) 



*o(/5) = M*o(/5) (4) 



where /i is the chemical potential and is chosen to be real 
and subject to the normalization condition 27rJp°°(ip/9'I'g — 1, 
so that no{p) ~ {N / L)'^'^{p) is the ground state density. The 
excited states are plane waves along z, with wavevector k. 
Different branches of such longitudinal waves exist, character- 
ized by the number of nodes in the radial direction, n. When 
the condensate is weakly excited in one of these modes, its 
order parameter can be written as 



*(p,z,<) 



(5) 



where 



+ <fe(p)e-''(^^— *)] . (6) 

Inserting this expression into Eq. Q one gets the Bogoliubov 
equation 



A. Bogoliubov excitations in trap 

Here we consider a condensate that is unbound along z 
and radially confined in the harmonic potential V{p) = 
(l/2)TO(jjp/9^. In this geometry one can better distinguish the 
different roles played by radial and axial degrees of freedom. 
The order parameter of the ground state only depends on p 



+ 9nn{p) \ I Unk \ _ _ I Unk 



where 



TuOnk 



Hp^-^ + Vip) + 2gnoip)-p 



Vnk 



(7) 



(8) 
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The quasiparticle amplitudes Unk and Vnk obey the following 
orthogonality and symmetry relations 

/ dp 2TTp{UnkU*n'k' - VukV^.k') = S„n'S{k - k') (9) 
/ dp p{UnkVn'k' - VnkUn'k') = . (10) 

Equations (0} and (|7j can be solved numerically to get ^fo, 
Unk, Vnk and LUnk 1 12, 24, 25, 26]. 

In Fig.|6lwe show the spectrum of axial phonons in a cylin- 
drical condensate which simulates the behavior of the elon- 
gated condensate discussed in section II. The relevant param- 
eter that characterizes the solutions of the above equations is 
77 = pTp/ifii^p) — [4aiV/i]^/^, where /ixF is the chemical 
potential in the Thomas-Fermi limit. The condensate in Fig.|6l 



has 77 = 9.1. Similar spectra were already discussed in details 

in Ref. LlZl . The form of the quasiparticle amplitudes Unk{p) 

and Vnk (p) for some states on the first two branches is shown 

in Figs. and |8](i = plots, first column). The modes on 

the lowest branch have no nodes in the radial direction; those 

in the first branch have one node, and so on. One also notes 

that the excited states have a nonvanishing amplitude in the 

center of the condensate only at low fc, while at higher k they 

are mainly located in the low density region near the surface 

iIhIIzt}! . Finally, we observe that, as expected, the amplitude 

V is of the same order of u in the phononic regime at small k, 

while it is negligible when k is larger than In our case 
^-1^^1/2 1^3 -1. 




FIG. 7: Time evolution of the functions \u{p) \ (solid lines) and \v{p) \ (dashed lines) obtained from Eq. <18> . At t = these functions coincide 
with the Bogoliubov amplitudes itofc and uofc , solutions of Eq. for the modes (a), (b), and (c) in the spectrum of Fig.|6| The radial coordinate 
p is given in units of ap. The Thomas-Fermi radius of the condensate is i? = 4.27ap. 



B. Free expansion and scaling 

Now let us suppose that the trapping potential is switched- 
off at i = 0. One can introduce a rescaled order parameter \^ 
asl28J 



where T is defined as r(t) = dt' / {t' ) . The scaling pa- 
rameter b{t) obeys the equation 



(12) 



imp'^b{t) iprit) 
2hb{t) h~ 



(11) 



with 6(0) = 1 and 6(0) = 0, whose solution is 



b{t) = [1 + u^t 



2 .21I/2 



(13) 
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With this choice, the GP equation ([Q for the rescaled order 
parameter ^{p, z, t) becomes 

f;2 n2 



A major property of this equation is that, if the rescaled or- 
der parameter at i = coincides with the t- and z-independent 
solution of the stationary GP equation (|4}, then it remains t- 
and z-independent solution of the same equation forever. Let 
us call 5*0 (p) this stationary solution. The time evolution of 



true order parameter is entirely accounted for by the scaling 
parameter b{t) given in ( I13t . The stationary "folp) thus corre- 
sponds to an expanding order parameter with constant shape 
in the rescaled co-ordinate p/b{t) and with a density that de- 
creases in time as n(p, i) {N/L)-^l{p/b{t))/b'^{t). This 
behavior is a particular case of a more general scaling prop- 
erty of the GP equation, which is exact in a two-dimensional 
geometry and is valid for a generic time-dependent trapping 
frequency 0Jp{t) |28, 29, 30]. Here we want to apply Eq. MA\ 
to describe also an expanding condensate whose initial con- 
figuration, at i = 0, is a superposition of the ground state and 
some excited states. 
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FIG. 8: Same as in Fig.0 but for the modes (d), (e) and (f). At t = the functions \u\ and \v\ coincide with the Bogoliubov amplitudes uik 
and vik, solutions of Eq. with one radial node. The position of the t — Q node is shown as a black dot on the horizontal axis. 



C. Rescaled Bogoliubov-like equations 

Let us take the rescaled order parameter in the form 
^{p, z, t) = *o(p) + ^*(p, z, t), where >Po is the (real) solu- 
tion of the stationary GP equation and no(p) — {N / L)^q{p) 
is the rescaled density. The linearized version of Eq. (I14t is 



ihdt5^{p, z, t) 



+ 



b^{t) 2m 



(5*(p,z,t) 
5-^*{p,z,t), (15) 



with 

Hp = -—^ + V{p) + 2gho{p)-p. (16) 
2m 

Now let us write (5^ in the form 

S^ip,z,t) ^ L-^/^[a{p,t)e''''' +v*{p,t)e-''''] . (17) 

Inserting this expression into Eq. il5i one gets the following 
equations for new amplitudes u and v 



ihdf 



b^t) " 



(18) 
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where 

These equations describe the behavior of small deviations 
from the stationary rescaled order parameter 5*0 . At t < 0, 
when b{t) ~ 1 and ^I* = ^l*, the t- and p-dependence of u and 
V can be factorized as in Eq. (|6j and the rescaled Bogoliubov- 
like equation ( I18> reduces to the Bogoliubov equation for 
the excitations of the trapped condensate. One can select one 
of these excitations, with given k and n, as initial values of u 
and V and then solve numerically Eq. (I18> for the free expan- 
sion at t > 0. Typical examples are given in Figs. and |8] 
where we show the time evolution of \u{p) \ and \ v{p) \ for six 
different excited states. 

During the expansion u and v exhibit different behaviors: 
u decreases in time at low k, while it remains almost constant 
at large k; vice-versa, v always decreases, but it decreases 
faster at high k. Let us define the normalized radial aver- 
age = Jdp p\dip,t)\y Jdpp\vip,0)\\ InFig.|9lwe 

plot (I^P) as a function of t for different k (solid lines). At 
large fc, all curve approach an universal behavior: v decreases 
in the typical timescale of the expansion, uj^^, which is the 
time needed for the mean-field energy to vanish. Conversely, 
at low k it decreases much more slowly, approaching a finite 
value for t ^ oo. This is also evident in Figs.0and|S] where 
the function \ v\ for the lowest value of k remains well visible 
also at t ~ 5 and t = lOujp^, when the mean-field inter- 
action is certainly negligible. In this regime, the functions u 
and V have completely lost their initial meaning of quasiparti- 
cle amplitudes of Bogoliubov modes; they are simply the two 
Fourier components of a density and velocity modulation in 
the expanding "ideal" gas. 

From the knowledge of u and v one can also write the 
rescaled density variation Sh, which is given by 

6n — ^o\u + v\ cos[fc2 + phase({t + v)] , (20) 

while {u — v) is related to the axial current density. The ini- 
tial phase can be chosen in such a way that the excitation at 
t = is just a sinusoidal wave of wavelength 27r/fc, which 
moves along z with phase velocity uink/k. Then one can 
look at the time evolution of nodal lines or crests. Typical re- 
sults are shown in Fig. [TO] The results at fc = 0.5a~^ can 
be compared with those obtained from the GP simulations 
in the elongated condensate for the same k, given in Fig. |5] 
The qualitative agreement between the two calculations is no- 
ticeable, especially if one keeps in mind that the results of 
Fig. [TO] completely ignore the inhomogeneity and finite axial 
size of the elongated condensate of Fig. |5] The k ~ ^a^^ 
case is also interesting, because it shows that the velocity of 
the fronts quickly approaches the phase velocity of free par- 
ticles hk/ {2m). Finally, from this results it is also clear that 
the observation of wavefronts, with their phase and amplitude, 
gives access to the quantity |t2 + v\, while the total momen- 
tum transferred is a measure of (|{tp — |12|. It is worth 
noticing that Eq. ( I18> conserves the total momentum, which 
is proportional to j dp p{\u\'^ — \v\'^). 
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FIG. 9: Normalized radial average of \v\'^ as a function of t, in units 
of LUp^, for different values of k, in units of a" ^. Solid lines are the 
results obtained from Eq. <18> for the infinite cylinder. Dashed lines 
are the results obtained from Eq. <23> or, equivalently, from Eq. <48> 
and correspond to the behavior of | {) p in a uniform gas whose density 
is equal to the radial average of the density of the expanding infinite 
cylinder. Dot-dashed lines correspond to the adiabatic following ap- 
proximation <45> . For larger values of k all curves converge to the 
analytic result <46> shown as a dotted line in the bottom-right quad- 
rant. 

D. Radial motion: scaling of nodal lines 

Now we consider in more details the time evolution of the 
radial shape of u and v. In the two opposite limits fc = 
and fc — > oo they obey an exact scaling law. Let us consider 
first the case fc = 0, i.e., purely radial excitations with an 
arbitrary number of radial nodes. The quantities u and v are 
z-independent and Eq. ( I18> becomes 

As initial condition att — 0, one can choose u{p, 0) — Uno{p) 
and v{p, 0) — w„o(p)^ where w„o and w„o are the eigenfunc- 
tions of the operator Cp with eigenvalues hujno^ i-e., the fc = 
solutions of the Bogoliubov eigenvalue problem 0. Then 
an exact solution of Eq. (12 1> is found in the form u{p, t) — 
a{t)uno{p) and v{p,t) = a{t)vno{p), with a{Q) — 1. One 
finds 

ib^t)dta{t) = LJnoc^it) , (22) 

whose solution is a{t) ~ exp[— icj„oT(t)] with r(i) — 
dt' /[I + iOpt'^] = ujp^ata.n{Ljpt). The phase of u and v 
has a nontrivial t-dependence, corresponding to a true oscil- 
lation of the type cxp{iijJnot) only for a short time interval, 
shorter than the typical expansion time oj'p'^, while for longer 
times the oscillations are frozen. However, the modulus |u| 
and \v\ remains stationary at all times and the radial nodes of 
these functions remain fixed. This means that the density vari- 
ation 6n has nodes that expand radially with the same scaling 
law of the whole condendate, governed by b{t). 
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FIG. 10: Wavefront position \'s. time for k 
2ap^, obtained from the maxima of Sn defined in Eq. <20t with u 
and V given by Eq. <18> . As in Fig. |5| the wavefront position is 
plotted for p = (solid), 0.4 (dot-dashed) and 0.87? (dashed), where 
R is the Thomas-Fermi radius of the condensate. The slope of the 
thin solid and dashed lines is the Bogoliubov sound velocity, c = 
[gn{0)/{2m)]^^'^ and the phase velocity of free particles, hk/{2m), 
respectively. The inset shows the shape of two nearest fronts at t = 
(dashed lines) and t — lOu}^^ (solid lines). Each front is plotted 
from p = 0top = i? = 4.27ap. The distance between the two 
fronts at t = 0, in units of ap, is 2n/k. 



The same scaling is found for k ^ oo when, since v is 
vanishingly small at all times, Eq. (I18> reduces to an equation 
for the amplitude u only: ihb^{t)dtu = [Hp + h^k"^ / {2m)]u. 
If M at t = is an eigenfunction Unk of [Hp + fi^k^/ (2to)], 
then it scales as u{p, t) — exp[— ici;„fer(t)]u„fc(/9), as before. 

In the general case of a finite and nonzero value of fc, the ra- 
dial and axial degrees of freedom in Eq. ( I18> are coupled and 
the functions [u{p)[ and [v{p)[ are no more stationary. How- 
ever, as one can see in Figs. Q and |8] their time evolution is 
very slow, at least for the Bogoliubov modes with few radial 
nodes or none, and the positions of maxima and minima are 



almost constant. This is particularly evident for the minimum 
of |u| in Fig.|S] which remains always very close to the posi- 
tion of the radial node at t = (solid circle on the axis). 



E. Axial motion: Thomas-Fermi approximation and local 
wavefront velocity 

Here we concentrate on the axial motion of quasiparticles. 
We first observe that when 77 ^ 1 the ground state density 
of the infinite cylinder is well reproduced by the Thomas- 
Fermi approximation gno{p) — /ixF — V{p), which coiTe- 
sponds to neglecting the quantum pressure — [?i^/(2m)]Vp 
in the stationary GP equation 0. This also implies that 
gnoip) — ptf ~ V{p) during the expansion. We also notice 
that the amplitudes u and v of axial phonons with a small num- 
ber of radial nodes are smooth functions of p at all times. This 
means that the term (2m)] Vp can safely be neglected also 
in Eq. (I18> . if we are not interested in the details of the ra- 
dial motion. With these approximations, expression (|8j gives 



Hp = gfiQ and Eq. jl8> becomes 



ihdf 



= H{t) 



where 



H{t) = 



Pit) 



1 1 

-1 -1 



2 1,2 



h^k 

2m 



1 

-1 



(23) 



(24) 



The same equation can be written for the quantities u + v and 
u — V L33.I : 



ih,dt{u + v) 
iTidt{u — v) 



2 1,2 



f 2gha 
Wit) 



{u — v) 



2m 



(25) 

(u + v) . (26) 



These equations can be easily decoupled to get the equation 
of motion for u + v: 



dt{u + v) = ~il\t){u + v) 



with 



fc2 f2gna Ti'^k'^ 



2m \b^{t) 



2m 



(27) 



(28) 



The quantity is the frequency of Bogoliubov excitations in 
a uniform gas having time-dependent density ho/b'^{t). Re- 
calling definition ( I20t . one can identify the quantity 



n 

'=k 



gnp 
mb^{t) 



hk 

2m 



1/2 



(29) 



with the axial phase velocity of density modulations associ- 
ated with the expanding phonons. This velocity depends on 
p through hf){p). One can insert the Thomas-Fermi density 
profile of the infinite cylinder and the expression of b{t) to get 
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analytic results for the wavefront position v.?. time. Typical re- 
sults of this local density approximation are shown in Fig. II II 
and can be compared with those obtained from the numerical 
integration of the rescaled Bogoliubov equations, in Fig. ^| 
There is a good qualitative agreement both at low and high 
k, even though, as expected, the local density approximation 
overestimates the p-dependence of c, especially for the mo- 
tion near the surface of the condensate, where Thomas-Fermi 
approximation fails. 




FIG. 1 1 : Same as in Fig. llOl but using the local velocity <29> to get 
the wavefront position at each p and t. 



IV. ADIABATIC QUASIPARTICLE-TO-PARTICLE 
EVAPORATION 

In Eqs. ( I23t - (l24li the quantities no, u and v depend on p, 
but one can extract useful results also by averaging out this 
weak p-dependence. A simple approach consists of taking 
a constant rescaled density uq equal to the radial average 
(r7-o(p))p = no(0)/2- Then, one can choose the initial u and 
V as the p-independent quasiparticle amplitudes of phonons in 
a uniform gas of that density and for a given k. Finally, one 



can solve Eq. J23t for the time evolution of the rescaled am- 
plitudes. The resulting behavior of the function nor- 
malized to its value at f = 0, is shown in Fig. |9] for different 
k (dashed lines). The remarkable agreement with the results 
for (Iwp) obtained from the rescaled Bogoliubov equations 
( I18> tells us that the evolution of a quasiparticle in the inho- 
mogeneous cylindrical condensate is very similar to the one 
in a uniform gas having a density that decreases in time as 

An interesting reformulation of the same problem is found 
by considering Eq. ( I23> as the evolution of a two level system 
governed by a nonhermitian time-dependent Hamiltonian. Let 
us rewrite Eq. i23\ in the form 



ihdt\S^) = H{t)\S^) 



with 



(30) 



(31) 



Then let us recall definition j28> and introduce the quantity 
Q{t) such that TiilsinhO — gho/b'^{t) and ?ir2cosh6 = 
gno/b^it) + h^k^/{2m). Thus 



tanh0(t) = gfiQ gfiQ 



2 1,2 



2m 



and Hamiltonian (I24> can be rewritten as 

H{t) = hn{t) 



coshe(t) 
-sinhe(i) 



sinhe(i) 
coshe(t) 



(32) 



(33) 



This Hamiltonian is nonhermitian. At any given time t, 
it admits two real eigenvalues ±0 and the corresponding 
biorthonormal set of eigenvectors {|±)r, |±);}, where the 
"right" and "left" eigenvectors are defined by 



H{t)\±)r = ±hn{t)\±)r 

HHt)\±)i = ±nnit)\±)i 



A simple calculation yields 



^ J_ / V cosh e + 1 
~ V2\ -\/coshe-l 

_ J_ f -Vcoshe - 1 

~ V2 V Vcosh e + 1 
1 



V2 



%/ cosh + 1 , V cosh 0—1 



1 



;(-| = ( Vcoshe- 1, Vcoshe + 1 ) 

These vectors obey the orthogonality relations 

/(+l+>r = = 1 

/(+h). - ;(-!+),■ = 0. 

They also satisfy the relations 

l{±\dt\±)r = 

i{T\dt\±)r = -(i/2)ate(t) 



(34) 
(35) 



(36) 
(37) 
(38) 
(39) 



(40) 
(41) 



(42) 
(43) 
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with 



Notice that the u and v components of the eigenvectors 
(|36}-(|39} are such that \u\'^ - [wp = ±1 for |±) states, 
with eigenfrequency ±f7. At i = the state \+)r, cor- 
responds to a quasiparticle of wavevector k and frequency 
n{0) = {fcV(2m)[2gno + h^k^ / {2m)]Y^^, whose compo- 
nents, u = Uk and ii — Vk, are the usual amplitudes of a 
Bogoliubov mode in a uniform gas. Alt ^ oo the same state 
corresponds to a free atom with the same k and frequency 
i7(oo) = hk^ /(2m). The time evolution of the state \+) r is 
simply obtained by using the expression ( I36> and the defini- 
tions of and VL{t). For example, one gets 



gno/b^{t) + rfcV(2m) 

2hn{t) 



(45) 



which has the two correct limits -5^(0) = = [gno + 
h^P/{2m)]/{2hn{0)) - 1/2 01 and ^^(oo) = 0. In the 
large k limit, when h k^ /{2m) ^ gnp, Eq. ( I45t gives 



mgnp 

h^k%^{t) 



(46) 



Now, let us use the basis of instantaneous eigenvectors 
{|±)r} to project the state |(54') as 



the term gfi^) in VI is larger than Ti'k^ / (2m) at t 
remains larger also att ~ ' ~^ 



0, but it 



iOp when fe^ — 2. The condition 



(44) for adiabaticity thus becomes fc ^ 77 
?/ = 9.1 and the condition is fc » 0.3a„ 



-1/2, 



In our case 




FIG. 12: Ratio |c_ j^/|c+ 1^ at f ^ 00 obtained by solving Eq. <48l 
for different values of fc, in units oi aZ^ , with the condition c+ = 1 



and c- 



at t = 0. 



15*) ^c+|+>,, + c_h). 



(47) 



and assume that the initial state is just a quasiparticle of mo- 
mentum hk, that is c+ = 1 and c_ = at < = 0. The evo- 
lution of |(5\I>) is governed by the Schrodinger equation j30> . 
By using the new basis together with the properties (I40ll-(I44>. 
the same equation can be rewritten in the form 



idt 



c+{t) 



n{t) 

^dtOit) 



^dte(t) 
-n(t) 



c+(t) 

c-{t) 



(48) 

The numerical integration of this equation gives, of course, 
the same results of Eq. ( I30> (dashed lines for |wp in Fig.|9j 
for instance), but the use of a different basis makes the evap- 
oration mechanism more transparent. In fact, it allows one 
to point out the interesting situation that occurs when \5'^) 
remains close to \+)r at all times. This adiabatic following 
happens when the diabatic coupling between the two eigen- 
states is small, i.e., when the off-diagonal term \dtQ\ is much 
smaller than the frequency splitting 20, |,35J. Now, the func- 
tion |9f0(i)| vanishes at t = and t — > 00 and is always 
smaller than dth^ / (2b^) which has maximum value u!p/2 for 
t — LUp^. Conversely, the function n(t) decreases monotoni- 
cally from i7(0) to fik'^/ (2m). Thus the evolution is certainly 
adiabatic if fik"^ /(2m) ^ LOp/2, that means k ^ a/,^. This 
condition is somewhat too restrictive. In fact, the two func- 
tions 19(91 and 20 might be comparable only around t = uj^^ 
and the adiabatic evolution is ensured if 51 ^ at that time. 
Now, for k of the order of aj^ 



One can check whether the quasiparticle-to-particle conver- 
sion is adiabatic or diabatic by directly comparing the evolu- 
tion of the state with that of the instantaneous eigenstate 
\+)r. The difference is a measure of non-adiabaticity. For the 
quantity in Fig.|9]the comparison has to be done between 
dashed lines, that come from the evolution of \5'^) through 
Eq. ( I48> . and dot-dashed lines, that come from the evolution 
of 1+),- through Eq. ( 145 l l. As one can see, for k ^ O.Sa^ ^ 
the adiabatic following approximation is indeed an accurate 
description of the evaporation process. The dashed and dot- 
dashed curves are already indistinguishable for A: — 2a~^ and, 
for larger k, all curves collapse on the dotted curve, which 
represents the asymptotic law ( I46t . For low k, conversely, 
the diabatic coupling between the |+),. and |— ),. eigenstates, 
originating from the time-dependence of the scaling parame- 
ter b(t), is no more negligible and the state \S'^) significantly 
differs from \ +)r both at short and long times. A nice view of 
this coupling can be seen in Fig.^] where we plot the value of 
|c_ p/|c+p att 00 as a function of k. In the adiabatic case 
this ratio vanishes. At low k, conversely, it tends to 1. From 
the definitions ( I36t -( l37t . one can easily see that the quantity 



at < = 00 is equal to 



and corresponds 



and rj = ghQ/(hujp) > 1, 



to the ratio between the number of atoms moving with oppo- 
site momenta ±:Tik att ^ 00, as a result of the evaporation 
of an initial quasiparticle with momentum -\-hk. It should not 
be confused, vice-versa, with the ratio |t)fep/|itfep at t = 0, 
which is different and is related to ±k components the mo- 
mentum distribution of the initial in-trap quasiparticle 1 3^^. 
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V. GP SIMULATIONS REVISITED 

We are now ready to use the results of sections III and IV 
to interpret some relevant features of GP simulations in elon- 
gated condensates. 

Axial motion of wavefronts 

Differently from the infinite cylinder, the elongated conden- 
sate expands also in the z-direction. Let us consider Eq. ( I29> 
and apply a local density approximation in this form: 



hk 



gnojp, z) 
mbl{t)h,{t) \2mh,{t) 



-1 1/2 



(49) 

The slow expansion along z is here included through the scal- 
ing parameter hz. The new scaling laws are l28lE9ll36ll 



K{t) = 



(50) 



Thus the first term in c is the local drift velocity due to the 
axial motion of the "background" in which the excitations 
move. The second term is now p- and z-dependent through 
the rescaled density no(p, z). This simple model allows us to 
plot the wavefront position and speed for the elongated con- 
densate, by using the unperturbed GP density profile at t = 
as input and solving the two coupled equations for and hp. 
The position of the wavefronts, corresponding to the simula- 
tion in Fig.|3j are shown as dashed lines in the same figure. In 
practice, we take the crests of 5n si t — from the GP simu- 
lation and let them move with velocity c during the expansion 
for t > Q. The dashed lines are the same crests alt ~ 25^^^^. 
They nicely agree with the results of GP simulations (solid 
lines), except near the surface where the Thomas-Fermi ap- 
proximation is inadequate. 

One can also estimate the position in time of the crests of 
the linear density (upper part of Fig|3}- A simple way consists 
in propagating each front with a radially averaged velocity, 
which is obtained by replacing no(p, z) in Eq. j49> with the 
average {nQ{p, z)) p. The results are shown as solid lines in 
Fig. El where they are compared with the ones of the GP 
simulation (empty circles). The thin dashed lines represents 
the Thomas-Fermi axial size of the expanding condensate. 
The condensate has maximum density at z = 0. All crests 
move initially with almost the same positive velocity along 
z. Within a short time interval, when the mean-field is still 
active, those crests that move downhill (positive z) tend to 
be accelerated by the expanding background, with respect to 
their motion in an infinite cylinder. Those that move uphill 
(negative z) seem to be decelerated, and they can also invert 
their motion. The agreement with the full GP simulation is 
remarkable. This means that Eq. i49i accounts for most of the 
physics involved in the axial motion of density modulations in 
the expansion of low k phonons. 

Threshold for the appearance of a released-phonon cloud 
Looking at Figs. 1101 andll3lone understands why at low k 
the excited atoms do not form a separate cloud: they have not 
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FIG. 13: Wavefront position vs. time for the elongated condensate 
of Figs. 1 and 3. Points correspond to the crests of the hnear density 
obtained from the GP simulation. Solid lines are the results of the 
local density model, with local phase velocity given by |49j. Thin 
dashed lines represent the axial size of the expanding condensate. 
Length and time are in units of ap and LOp^, respectively. 



enough speed to reach the expanding boundary along z. In 
classical fluids this situation corresponds to the case of sound 
waves propagating on a background hydrodynamic flow. A 
boundary that moves faster than the speed of sound acts as 
an event horizon in a sonic black hole |37]. A possible re- 
alization of a stable black hole for long wavelength phonons 
in ring-shaped Bose-Einstein condensates has been recently 
proposed in Ref. 138]. Our configuration is nonstationary and 
the excitations are sound-like waves only at the very begin- 
ning of the expansion, in a time interval of the order of uip"^, 
when they are still well inside the condensate. Then, the wave- 
fronts reach an almost constant phase velocity, of the order of 
fik/ (2m) and the whole wave packet moves at group velocity 
~ fik jm. Thus, if the latter is smaller than the velocity of the 
boundary (upper dashed line in Fig.ll3>. most of the atoms re- 
main within the expanding condensate at all times. In the limit 
A ^ 1, Eqs. ( 15 Oi l yield the expression (7r/2)wpA^ZTF for 
the asymptotic velocity of the boundary, where Ztf is its ini- 
tial position in Thomas-Fermi approximation |29]. Thus the 
condition for the atoms to exit is hk/m > {tt / 2)uj p\^ Z^y ■ 
Since (l/2)mA-^w^Z|,p = /itf and 77 = fiTF /{^ojp), the 
same condition can be rewritten as fcop > 7rA(77/2)^/^. With 



the parameters of the condensate in section II, this threshold 
is around k ~ O.TGa^ ^, in agreement with the transition from 
the low k scenario of Fig. ^ to the high k scenario of Fig.|2 
The transition is not sharp, however, since the excited atoms 
can reach the boundary at different times and with different 
velocities so that, for a given expansion time, one can ob- 
serve a released-phonon cloud only partially separated from 
the condensate. 

Adiabatic and diabatic evaporation 

For the elongated condensate of the simulations in section 
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II one has 7y-i/2 < 7rA(?7/2)i/2 < r\^l'^ . The first value is 
the threshold for the adiabatic quasiparticle-to-particle evapo- 
ration discussed in section IV (see Fig. ll2> . The second is the 
threshold for the appearance of a separate released-phonon 
cloud. The third is the inverse of the healing length. This 
implies that when the initial excitations are single-particle 
modes, i.e., with fcop > r]^!'^ , the excited atoms always move 
out of the condensate. Instead, when the initial excitations 
are phonon-like modes, i.e., with fcop < ?7^/^, they can ei- 
ther move out or remain inside. However, if they move out as 
a separate cloud, the evaporation is certainly adiabatic. This 
means that a single quasiparticle of positive momentum fik 
gives rise to a single atom moving with the same positive mo- 
mentum, even though the initial quasiparticle state is made 
of correlated +fc and — fc components in the momentum dis- 
tribution of the atoms. It is worth noticing that, with appro- 
priate choices of A and ry, such that 7rA(r//2)^/^ < rf^l"^, 
one can have diabatic evaporation together with the appear- 
ance of released-phonon clouds. In this case, one expects to 
have atoms moving with both positive and negative momen- 
tum ±?ifc out of the condensate. 
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FIG. 14: Radial density profiles, in arbitrary units, of the released- 
phonon cloud at z = hkt/m, with t = 25u)p^, k — 2.5ap^ and 
three different values of lo. Points with dashed lines are the results of 
GP simulations with Vb = O.Shcup and ts = icop^ . Solid lines are 
the functions \unk + v^k] for an infinite cylinder with 77 = 9.1, the 
same k and n = (bottom), 1 (mid) and 2 (top). 



to po/b{t), but translated along z as a consequence of the ax- 
ial motion. If this is true, the shell structure in Figs.|21and|3 
is an indication that, even though we used a Bragg frequency 
in resonance with the lowest phonon branch with no nodes, 
some quasiparticles in the next branches, with one and two 
radial nodes, were also excited. In order to be more selective, 
one has to use Bragg pulses with longer duration ts fTl','T^. 
We thus repeat the same type of simulations of Figs.|2]and|4] 
for k = 2.5a~^, but with is = 4:Ujp^. The Bragg frequency 
is varied in such a way to excite phonons with 0, 1 and 2 ra- 
dial nodes. We find a significant change in the shell structure 
of the released-phonon cloud. In Fig. we plot a cut of its 
density distribution at z — hkt/m, with t — 25ujp^, which 
corresponds roughly to the position where the moving cloud 
has its maximum radial extension. The three density profiles 
obtained from GP simulations (points with dashed lines) are 
compared with the shape of the functions \unk + Vnk\ (solid 
lines), calculated at p/b{t), for the Bogoliubov modes of an 
infinite cylinder, with the same k and with n = 0, 1 and 2. 
The good agreement confirms that the nodal lines scale al- 
most exactly and that the observation of radial structures in 
the released-phonon cloud can be an efficient tool for the iden- 
tification of the initial in-trap quasiparticles. 

VI. CONCLUSIONS 

We have shown that the evaporation of phonons in a freely 
expanding condensate is an interesting quantum process with 
a rich variety of scenarios. Our analysis is based on the GP 
theory at zero temperature and in the linear response regime, 
where the concept of quasiparticle is well defined. We have 
thus neglected possible effects of thermal and quantum fluc- 
tuations, as well as those of nonlinear dynamics. This is justi- 
fied if one assumes that the population of the initially excited 
phonon-like mode is large enough to neglect fluctuations and 
small enough to remain in the linear regime. For condensates 
with 10^ to 10^ atoms or more, this condition is ensured by 
choosing the Bragg intensity and duration such to excite about 
1 to 10% of the atoms as in typical experiment. In this paper, 
we have provided explicit calculations for an elongated con- 
densate similar to that of current experiments, but the results 
obtained for the infinite cylinder are rather general. The ex- 
tension of this analysis to different sets of parameters can lead 
to different regimes that are also worth exploring. 



Radial scaling of nodal lines 

In section II, we have seen that the released-phonon cloud 
exhibits a "shell" structure (see Figs.|2land|4j. It can be rea- 
sonably assumed that the scaling of the radial motion, which 
is almost exact for the infinite cylinder, is still valid also for 
the elongated condensate. This means that an excited state 
having a radial node at a point po{z) at i = should produce 
a cloud of excited atoms having a minumum of density close 
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